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A Monte Carlo method to sample the classical configurational canonical ensemble is introduced. In 
contrast to the Metropolis algorithm, where trial moves can be rejected, in this approach collisions 
take place. The implementation is event-driven, i.e., at scheduled times the collisions occur. A 
unique feature of the new method is that smooth potentials (instead of only step-wise changing 
ones) can be used. Besides an event-driven approach where all particles move simultaneously, we 
also introduce a straight event-chain implementation. As proof-of-principle a system of Lennard- 
Jones particles is simulated. 
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I. INTRODUCTION 

The most commonly used methods for simulating par- 
ticle systems in accordance to classical statistical me- 
chanics are molecular dynamics (MD) and Monte-Carlo 
methods (MC) based on the Metropolis scheme [THS]. 
For systems, such as hard-sphere systems, with impul- 
sive interactions a time-driven MD approach does not 
work and an event-driven approach can be used. In fact, 
the pioneering work of Alder and Wainwright used an 
event-driven molecular dynamics (ED-MD) scheme [5J. 

In the Metropolis MC scheme trial moves are either 
accepted or rejected. In highly concentrated systems the 
acceptance rate can be very low and simulating using 
MD requires very small time-steps. In dilute systems 
the time-scale in MD or step-size in MC is determined 
by the molecular collision process and simulation time 
is wastefully spend on flying through empty space. In 
both cases an event-driven approach can speed up the 
computation. 

ED-MD can be generalized to hard-spheres to poten- 
tials build up by a sequence of steps [7] ■ Clearly in this 
case an event takes place at each step. The method we de- 
rive in this paper differs in several aspects from ED-MD: 
Collision-events are determined by means of a stochastic 
process. Potentials are not necessarily step-wise. There 
is no exchange of kinetic and potential energy. In fact mo- 
mentum is not relevant and the configurational canonical 
is sampled directly. 

Instead of rejecting moves as in the Metropolis scheme 
a collision takes place. There is quite some freedom to 
model a collision event. One possibility is to model it as a 
Newtonian collision. Another possibility is to move one 
particle at a time where, at collision, another particle 
takes over. This is similar to the straight event-chain 
collision in hard-sphere simulation 19] . 

On an algorithmic level there is some similarity with 
kinetic (or dynamic) MC [TU] and n-fold way MC simula- 
tions jTTJ [T2] . In these methods there is a finite number 
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of (classes of) moves modeled as Poisson processes. Us- 
ing the rates corresponding the these Poisson processes 
the moment in time a next event occurs can be com- 
puted. The efficiency of these methods is determined by 
the fact that the number of moves is finite which is not 
the case in a particle system. Also in these kinetic MC 
simulations nothing happens in between two subsequent 
events. In the method we will outline below, however, 
particles will move linearly in between subsequent colli- 
sion events. Therefore even when no events occur the 
system is evolving. The present method, which is sur- 
prisingly simple, is a unique event-driven Monte Carlo 
method. 



II. AN EVENT DRIVEN STOCHASTIC 
SCHEME 

The prototypical Monte-Carlo scheme for sampling a 
configurational canonical distribution generates "moves" 
from an old state, x" ld , to a new state, x" ew , according 
to a conditional probability density T(x™ ow |x™ ld ). The 
transitional probabilities are forced to obey the detailed- 
balance relation, 

T(y"|x") exp[-/3C/(x")] = T(x"|y") exp[-0 U(y n )}. 

In the Metropolis scheme we decompose the transition 
probability density as, 

T(y«|x") = acc(y",x n )a(y"|x«), (2) 

where a(y™|x ra ) is the probability density for generating a 
trial move from x™ to y™ and acc(y ra , x") the probability 
that this move is accepted. The Metropolis form for the 
acceptance probability equals, 

acc(x£ ew ,x£id) = mm(l,exp[-/?At/]) , (3) 

if a(y ra |x n ) = a(x"|y") Vx n , y". When a move is not 
accepted the positions remain unchanged: x" ow := x™ ld . 

Now let's consider a simple one-dimensional potential 
step of height AU. In this case detailed balance, Eq. ([I]), 
can be obeyed in a different way. Instead of rejecting a 
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FIG. 1. (Color online) A trial move that moves upward to 
the higher energy state may give rise to a collision. 



move, if a random number is below the Metropolis accep- 
tance probability, it will collide. So, let's consider a trial 
move from a position x id to i ncw . If both positions are 
at the same side of the barrier the move will be accepted. 
If the move descents the barrier, i.e., AU < then the 
move is also accepted. If the move is up the barrier, i.e., 
AU > it will only sometimes be accepted. If it is not 
accepted, it is not rejected but the path is changed by 
means of a collision against the "wall" of the barrier (see 
Fig. I]) . Clearly a position ir now that is on the other side of 
the barrier as a; id can only be sampled if no collision has 
taken place. For the probability that no collision occurs 
we USe Eq. Pno-coll^newj^old) = acc(£„cw,2:oid)- 

Next let's consider a number of potential steps in a 
sequence. For a trial move from x id to x ncw we compute 
the probability to not collide at each individual barrier 
that is crossed by means of Eq. ^. In this case the 
probability to still have not experienced any collisions 
when reaching position x new equals, 



min (1, exp[— p AUi]) 

i 

= exp \-P max( AC/i , 0) 



(4) 



where the index i labels the barriers crossed when mov- 
ing from x id to Xncw For every change in potential we 
decide to count it or not depending on the fact whether 
it is increasing the potential energy or not. Going down 
the barrier is free, every uphill motion counts and accu- 
mulates until a collision becomes inevitable (or until the 
potential does not grow anymore). 

We could approximate a continuous potential by a se- 
quence of barriers and do our calculation accordingly but 
we will proceed differently. If we take the limit to indef- 
initely small potential steps we obtain, 
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FIG. 2. (Color online) The upper-left graph shows the rele- 
vant part of the potential for motion to the right starting at 
x = —1.5 and the upper-right graph for motion to the left. 
The cumulative probabilities to be not collide are shown in the 
second row of graphs. In the lower-left corner a typical time- 
series is depicted. The symbols indicate equidistantly spaced 
points along the s axis. These points sample the canonical en- 
semble as shown in the lower-right graph (solid line). When 
making a histogram of the collision points one finds the dashed 
curve. 



U a , where the subscript a is just a label to identify the 
potential (which is useful for reasons that will become 
apparent below). 

In practice computation of the integral is trivial if one 
has an expression for the potential t/ Q (x"(s)). One needs 
to know the location maxima and minima of the poten- 
tial along the path, x n (s), to be able to extract increasing 
contributions only. With this accumulative probability 
the position at which the particle does collide can be de- 
termined as follows: Draw a uniform number, u, between 
and 1. The collision takes place at the time, s, for which 



Dii(s), or equivalently, 
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maxf ^-U a 
\ds 



™(s)),0) ds = -kTlnu. 



(6) 



,(s) = exp 



lliax(^t/ a (x n (g)) > 0) lis 



ds 



(5) 

which is the conditional probability that a particle mov- 
ing in a linear motion from x n (so) to x ra (s) did not expe- 
rience a collisions along the way. Here we presented the 
formula for a general n-particle system, and a potential 



A. 1-D proof-of-principle 



To prove that the scheme correctly works in practice we 
consider the motion in a harmonic well: U = \ x 2 . Here 
we use dimensionless units kT = 1 and the characteristic 
length scale equals 1. 
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FIG. 3. (Color online) A particle moving in two potentials 
indicated by the two sets of equipotential contours. 



III. A 3-D MULTI-PARTICLE SYSTEM 

Let's consider a particle system. Here the total poten- 
tial can decomposed as a sum of potentials: V — J2 a U a - 
In Fig. [3] a particle moving in the fields of two potentials 
is shown. Now, assume for a moment that the potentials 
do not increase smoothly, but stepwise at every depicted 
equipotential contour. Using the same reasoning as be- 
fore, at every step that is crossed by the path of the 
particle a collision can take place. The probability that 
a path of the particle crosses a step of both potentials 
exactly at the same time is zero. Therefore, in the step- 
wise case, it is clear that the influence of each potential 
XJ a can be considered separately and this remains valid 
in the limit of smooth potentials. For each potential V a 
individually Eq. ^ can be used. 



A. Collision rules 



The motion of x is linear dx/ds — v (constant v) and 
at a collision: v := —v. Note that after the collision also 
dV/ds has changed sign. At this point, say at position 
afcou and time s co ii, we proceed with the linear motion 
and determine the new cumulative probability not to col- 
lide by means of Eq. ^ and integrating from an initial 
position Xcoii. Using this new cumulative probability the 
next collision is determined by means of solving Eq. (JTjj) 
with s = Scoii- 

To illustrate the process, in Fig. [2] the particle starts 
to move at x = —1.5. First, the collision "time" and po- 
sition are determined, then the collision is performed by 
reversing the "velocity". Here quotation marks are used 
because not "time", but contour length s is the relevant 
parameter. The "velocity" does not have physical signif- 
icance, e.g., as used for a kinetic energy. It is, however, 
more intuitive to speak in terms of time as the variable 
that parameterizes the path. 

To generate the canonical ensemble the positions, x, 
need to be sampled at equidistant points in time. If we 
define a time-step, say As, at every time s n — nAs the 
distribution is sampled. In the lower-left graph of Fig. [2] 
the data points corresponding to As = 2 are shown in 
the time-series. When collecting these points to form a 
histogram the correct canonical ensemble is sampled as 
is shown in the lower-right graph. It is a rigorously valid 
procedure, obeying detailed balance, if a new velocity v 
is drawn from a probability distribution, which is even 
in v, at equidistantly spaced times. In the series gener- 
ated to produce the bottom graphs, however, we do not 
do this and just proceed along the path until the next 
collision occurs. The dynamics has enough inherent ran- 
domization to cause ergodicity. The velocities are -1 or 1 
with equal statistical weight and clearly not distributed 
according to, e.g., a Maxwell-Boltzmann distribution. 



Let the particles in the system move with constant ve- 
locity, v n = (vi, v 2 , . .., v n ). If a collision due to poten- 
tial V a takes place at time s co n the velocity after collision 
changes as, 

v"(4on) = (I-2P Q )-v"(s- oU ), (7) 

where P a is a projection matrix (P a • P Q = P a ). A 
general form for the projection operator is, 

_ M-W a W a 
a VV a -M-VV a [ > 

Here the potential gradient indicates the direction normal 
to the equipotential surface of V a . One can verify that 
the collision leaves the scalar v™ • M _1 ■ v" invariant. 

A possible simulation protocol proceeds as follows: 
Draw velocities from a Gaussian distribution with a co- 
variance matrix proportional to M. Next run the event- 
driven collision scheme for a time-interval As = 1. Lastly 
redraw the velocities and repeat. This scheme gives rise 
to a Markov chain that obeys detailed balance. In our 
simulation results we find that, in fact, the velocities do 
not need to be redrawn. In appendix ^ we provide a 
proof that the algorithm indeed samples the configura- 
tional canonical ensemble. 

In the case that a pair-potential acts between particles 
1 and 2, J7 a (x n ) = J7(xi,X2), the potential gradient row- 
vector only has non-zero entries for particles 1 and 2. 
For the simulations we made the simple choice M = I. 
For pairwise central potentials, V a (x n ) — V(\x2 — xi|), 
we find P a ,ij = 5(^1^1 + <5i2^2 - 6aSj2 ~ fe^i) e r e r , 
with e r the radial direction vector e r = (x 2 — xi)/|x2 — 
Xi|. This is a formal notation equivalent to an elastic 
Newtonian collision between two particles of equal mass. 

The simulation protocol is very similar to event-driven 
MD [13]. Initially for all possible pairs a possible colli- 
sion event is computed and stored in a priority queue. If 
the collision that involves particles i and j pops up it is 
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handled. Now all previously computed collisions involv- 
ing i or j become invalidated and are removed from the 
queue. So, for all pairs i — k and j — k new collision times 
need to be computed similarly as in Eq. ^ by inverting 
Eq. ([5]). From a computational point of view it is most 
efficient to perform the updating asynchronously, i.e., the 
particles are moved only at the moment they participate 
in a collision, otherwise the positions remain fixed at the 
spot the last collision occurred. However, to generate the 
statistics we need to sample the system at equidistantly 
spaces time-intervals s stamp = nAs. We also schedule 
these time stamps, such that at every time s s tamp the 
positions of the particles, x™(s sta mp), can be computed. 

B. Straight event-chain collisions 

It has recently been shown that straight event-chain 
updates can be very efficient for concentrated hard core 
systems [HUH]- This makes implementation simpler than 
for the scheme outlined above because no event-queue is 
needed. Hence we tested also this scheme. 

If a particle, i, that moves with a velocity, v$ = v, col- 
lides with a particle, j, it stops (v$ := 0) and the other 
particle takes over (vj := v). The motion with colli- 
sions continue until As = 1. It was found that, when 
this scheme was performed non-reversibly, e.g., by giv- 
ing particles either one out of three possible velocities: 
v = ve x , v = ve y or v = ve z (and not the negative di- 
rection), the speed up was significant. The reason is that 
the dynamics is non- diffusive. Clearly in this case de- 
tailed balance is not obeyed but for hard sphere systems 
it was found that the correct configurational canonical 
ensemble is sampled. 

C. Lennard- Jones interaction 

As a second example we will have a look at the 
truncated-shifted Lennard- Jones interaction. 

((;)"-(;)*) 

r ( 9 ) 

fjtmnc^s \Uu{r)-U u (r c ), forr<r c 
LJ K ' jo, otherwise 

When particles move towards each other the pair- 
potential differences increases when it is in the repulsive 
regime. Once beyond the point of closest approach, and 
inside the repulsive regime, the motion is downhill and 
no collision can occur there. If particles move away from 
each other the potential difference increases if the parti- 
cles is inside the attractive regime of their pair-potential. 
Here a collision might occur. The parts that contribute 
to an increase in the collision probability as computed 
from Eq. ^ are illustrated in Fig. [4] 

Fig. [5] shows the radial distribution function (RDF) 
for the density p = 0.317 and T = 1.085 in LJ-units (and 




FIG. 4. (Color online) Particles interact with the central bead 
by means of a truncated-shifted Lennard- Jones interaction. 
The inner dash-dotted circle indicates the location of the po- 
tential minimum. The outer dot-dashed circle indicates the 
location of the cutoff radius. The bold solid pieces of the 
particle trajectories indicate the parts where the potential in- 
creases when the motion proceeds. These parts contribute in 
Eq. jBJ. In the other sections of the paths no collision can 
occur. 

1000 particles) for the case r c = 2.5a. This is the critical 
point of this truncated-shifted L J-potential [14] . We have 
chosen this point, instead of e.g. a liquid state point, 
because here there is a clear influence of the attractive 
part of the potential on the RDF. 

This RDF has been determined with the event-driven 
scheme where all particles move simultaneously for two 
cases. In the first case velocities are periodically redrawn 
from the correct Gaussian distribution. In the second one 
the velocities are never reset. We also implemented both 
a reversible and irreversible version of the straight event- 
chain method. As a check the RDF was also computed 
using the Metropolis scheme. All curves are identical 
within statistical errors. The maximal absolute deviation 
amongst the presented curves is 0.006 near r = 1.1. 



IV. DISCUSSION 

The event-driven rejection-free MC method outlined in 
this paper was successfully applied to a Lennard- Jones 
fluid. We only considered pair-potentials. If one wants 
to simulate a molecular system also angle and torsion 
potentials need to be considered. The collision rule such 
as defined by Eq. ^ can also be used for these kind of 
potentials. In that case 3 or 4 particles are involved in 
a collision, but solving Eq. ([6| will require some more 
computational effort. The generalization of the straight 
event chain collisions to these kinds of potentials seems 
less trivial. 

A priori it is not clear if for molecular systems the 
new method is less efficient than MD or not. As demon- 
strated by the harmonic well example Fiq. [2j the mo- 
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Upon collision the velocity changes according to 0. As 
a shorthand notation for the collision operator we will 
use R Q = (I — 2P Q ). Two relevant properties of this 
operator are 

R a • R a = I and R„ • VU a = -VU a . (A2) 

After a collision the velocities, v™, become R Q • v™ and, 
vise- versa, R Q • v™ changes into v™. 

The change of p(x n , v" , s) with time has three contri- 
butions: streaming, creation of states with velocities v ra 
due to collisions and annihilation of states with velocities 



FIG. 5. The solid line shows the RDF of the truncated-shifted 
Lennard- Jones potential with r c = 2.5 at p = 0.317 and T = 
1.085 (LJ-units). The curves computed with the rejection-free 
method and the Metropolis method are identical. 



tion goes from one side of the potential to the other. In 
MD one needs to resolve the oscillating motion by us- 
ing sufficiently small time-steps. The time that is won 
in this way can be spend on the more involved compu- 
tation of computing events and maintaining the event- 
queue. Although MD simulation tools are quite mature 
the algorithms for event-driven simulations are still being 
improved [151 116) . The method presented in this paper 
widens the realm of possible applications of the event- 
driven particle because a large class of potentials can be 
handled now. It remains to be seen if the application of 
the method is suited for niche applications only, or that 
it can rival with MD and Metropolis-MC for general pur- 
pose molecular simulations. 



Appendix A: Proof of correctness 

In this appendix we will proof the validity of the 
rejection-free scheme. We will do this by demonstrating 
that the canonical distribution is the invariant distribu- 
tion of the dynamics of the system. 

The probability distribution to have at time s a particle 
system with positions x™ and velocities v™ is denoted by 
p(x n ,v n ,s). The total potential of the system is given 
by ^2 a U a . The probability to have not collided with 
a potential U a is given by ([5]). The probability density 
per unit time to collide with potential U a when in state 
(x™,v™) equals 



d_ 

ds 



p(x n ,v",s) = -v n - Vp 



Pcoll , Pno— ' > 

as 



= /3max(v™ • VU a ,0 



u, tt (a) =/3max(^[/ Q (x"( S )),0) 

(Al) 



+ /? m ax(V ■ Rl • VU a , o) p(x n , R Q • v™, s) 

a. 

- /3^max(v" • VC/ Q ,o) p(x",v n ,s) (A3) 



In this equation the gradient operator denotes differen- 
tiation towards positions only and not towards veloci- 
ties. From the second relation in (A2) we find that 
v™ • Rj • VC/ a = — v" ■ 'VUa. Furthermore, from the 
definition of the projection operator, ([8]), one can derive 
that the scalar (v n ) ■ M _1 ■ v™ is an invariant of the colli- 
sion operator R Q for any a. Therefore, if we assume the 
form p(x n , R Q v™, s) = p x (x n , s) Hv n ■ MT 1 • v"), we find 
that for the collision terms of d A3|) 



f3^2 max(v" • Rj • VU a , OJ p{x n , R a ■ v™, s) 

a 

-/3^max(v n • VC/ Q ,o) p(x n ,v n ,s) 
= /3^(max(-v" • VU a , o) - max(v" • VU a , o))p a • / 

= ^(E v "- vu a ) Px -f. 

a 

(A4) 



Using this result we find that for a canonical distribution, 
p x = Z^ 1 exp|— /? U a ], the streaming part and the 
collision terms in ( |A3[ ) cancel. This concludes the proof 
that the configurational canonical ensemble is indeed an 
invariant distribution of the dynamics generated by the 
rejection-free method. 



[1] K. Binder and A. Baumgartner, Applications of the in current physics (Springer- Verlag, Berlin ; New York, 

Monte Carlo method in statistical physics, 2nd ed., Topics 



6 



1987) pp. xvi, 341 p. 

[2] M. P. Allen and D. J. Tildesley, Computer simulation 
of liquids, Oxford science publications (Clarendon Press 
; Oxford University Press, Oxford England New York, 
1989) pp. xix, 385 p. 

[3] D. Frenkel and B. Smit, Understanding molecular simu- 
lation : from algorithms to applications, 2nd ed., Compu- 
tational science series (Academic Press, San Diego, 2002) 
pp. xxii, 638 p. 

[4] D. C. Rapaport, The art of molecular dynamics simula- 
tion, 2nd ed. (Cambridge University Press, Cambridge ; 
New York, 2004) pp. xiii, 549 p. 

[5] W. Krauth, Statistical mechanics : algorithms and com- 
putations, Oxford master series in physics (Oxford Uni- 
versity Press, Oxford, 2006) pp. xii, 342 p. 

[6] B. J. Alder and T. E. Wainwright, J Comput Phys. 31, 
459 (1959). 



[7] G. A. Chapela, S. E. Martinezcasas, and J. Alejandro, 

Mol Phys. 53, 139 (1984). 
[8] E. P. Bernard, W. Krauth, and D. B. Wilson, Phys Rev 

E 80 (2009). 

[9] E. P. Bernard and W. Krauth, Phys Rev Lett. 107 
(2011). 

[10] K. A. Fichthorn and W. H. Weinberg, J Comput Phys. 

95, 1090 (1991). 
[11] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J Comput 

Phys. 17, 10 (1975). 
[12] M. L. Guerra, M. A. Novotny, H. Watanabe, and N. Ito, 

Phys Rev E 79 (2009). 
[13] D. C. Rapaport, Prog Theor Phys Suppl. 178, 5 (2009). 
[14] B. Smit, J Comput Phys. 96, 8639 (1992). 
[15] S. Miller and S. Luding, J Comput Phys. 193, 306 (2004). 
[16] M. N. Bannerman, R. Sargant, and L. Lue, J Comp 

Chem. 32, 3329 (2011). 



